clear; if ispc; dbstop if error; end; rng('default');

input_folder =  [fullfile(pwd, '..') filesep];

output_folder = input_folder;

%% load estimates

load([input_folder 'fc_covariates2016.mat'], 'dev_profit_one','dev_profit_all','dev_profit_one_boot','dev_profit_all_boot', ...
    'in_state_sim', 'craft_sim', 'entry_sim', 'pop_profit');

[N_mkt, N_dev, N_boot] = size(dev_profit_one_boot);


disp(['N_boot = ', num2str(N_boot)])


parpool('local', maxNumCompThreads);

%% setup

use_scope = false;
est_corr = false;  % whether to estimate the correlations
N_dev2 = 200;    % how many product pairs to use for estimation 

if ~use_scope && ~est_corr
    reps = 500;
else
    reps = 200;  % number of simulants to calculate the critical value
end

sig_level = 0.05;  % significance level
N_test = 20;  % number of tests to determine the minimum of the test statistics
N_grid = 200; % number of trial points
N_perturb = 50; % # of forward perturbations


est_years = 2016;

% market size cutoffs for defining small/median/large markets
pop_profit_cutoff = [0 1e5 5*1e5 Inf];

% profit cutoffs for defining ivs
profit_cutoff = [2.5 10]; % this is [30 120]/12

%% estimation
    
% construct covariates
fc_x = nan(N_mkt, N_dev, 5); % baseline fixed cost covairates

fc_x(:, :, 1) = in_state_sim; % whether the brewer is in state (CA)
fc_x(:, :, 2) = craft_sim; % whether the product is craft
fc_x(:, :, 3) = repmat(pop_profit<=pop_profit_cutoff(2), [1 N_dev]);% the market is small
fc_x(:, :, 4) = repmat(pop_profit>pop_profit_cutoff(2) ...% the market is medium
    & pop_profit<=pop_profit_cutoff(3), [1 N_dev]);
fc_x(:, :, 5) = repmat(pop_profit>pop_profit_cutoff(3) ...% the market is large
    & pop_profit<=pop_profit_cutoff(4), [1 N_dev]);

ind_mktsize = [4 5]; % index of fc_x that correspond to market dummies with a different variance parameters
ind_mktsize_scope = [4 5];

% construct instruments
fc_iv0 = nan(N_mkt, N_dev, 4);
fc_iv0(:, :, 1) = in_state_sim;
fc_iv0(:, :, 2) = craft_sim;
fc_iv0(:, :, 3) = repmat(pop_profit<pop_profit_cutoff(2), [1 N_dev]);
fc_iv0(:, :, 4) = repmat(pop_profit<pop_profit_cutoff(3), [1 N_dev]);

% reshape into product-market level obs
N_fc_x = size(fc_x, 3); % # covariates in the fixed cost function
N_fc_iv = size(fc_iv0, 3);% # instruments

fc_x_tr = reshape(permute(fc_x, [2 1 3]), [N_dev*N_mkt N_fc_x]);% fixed cost covariates
fc_iv0_tr = reshape(permute(fc_iv0, [2 1 3]), [N_dev*N_mkt N_fc_iv]);

entry_sim_tr = entry_sim'; entry_sim_tr = entry_sim_tr(:); %product entry: N_dev*N_mkt x 1. Each N_dev elements correspond to a market.    
dev_profit_one_tr = dev_profit_one'; dev_profit_one_tr = dev_profit_one_tr(:);% ub of change in variable profit : N_dev*N_mkt x 1
dev_profit_all_tr = dev_profit_all'; dev_profit_all_tr = dev_profit_all_tr(:);% lb of change in variable profit: N_dev*N_mkt x 1


dev_profit_one_boot = reshape(... 
    permute(dev_profit_one_boot, [2 1 3]), [N_mkt*N_dev N_boot]);


dev_profit_all_boot = reshape(...
    permute(dev_profit_all_boot, [2 1 3]), [N_mkt*N_dev N_boot]);


% indices of parameters
ind_para_fc_x = 1 : N_fc_x;
ind_para_sig = N_fc_x+1 : N_fc_x + length(ind_mktsize) + 1;% sigma parameters
para_ind_end = max(ind_para_sig);
    
if use_scope
    ind_para_scope = para_ind_end + 1 : para_ind_end + length(ind_mktsize_scope) + 1;
else
    ind_para_scope = [];
end

if est_corr
    ind_para_corr = para_ind_end + 1;
else
    ind_para_corr = [];
end

% define ivs    
g_ar_tr = gen_iv(dev_profit_all_tr, dev_profit_one_tr, profit_cutoff, fc_iv0_tr);

% add first stage noises
entry_sim_boot = repmat(entry_sim_tr, [1 N_boot]);
fc_x_boot = repmat(fc_x_tr, [1 1 N_boot]);
g_ar_boot = zeros(N_mkt*N_dev, size(g_ar_tr, 2), N_boot);
for nb = 1 : N_boot 
    g_ar_boot(:, :, nb) = gen_iv(dev_profit_all_boot(:, nb), ...
        dev_profit_one_boot(:, nb), profit_cutoff, fc_iv0_tr);
end

% generate additional covariates for estimating economies of scope
if use_scope
    scope_script;
else
    brewer_entry_tr = [];
    max_profit_top_vec = [];
    min_profit_top_vec = [];
    max_profit_top_vec_boot = [];
    min_profit_top_vec_boot = [];
    brewer_fc_x_tr = [];
    brewer_g_ar_tr = [];
    brewer_g_ar_tr_boot = []; 
    fc_err = [];
    N_brewer_dev = [];
end

if est_corr
    corr_script;
else
    entry_dev2 = [];
    dev2_ind_f1_tr = [];
    dev2_ind_f2_tr = [];
    corr_g_ar = [];
    corr_g_ar_boot = [];
    corr_sim = [];
end

clear entry_sim in_state_sim craft_sim pop_profit dev_profit_one dev_profit_all dev_profit_one_boot0 dev_profit_all_boot0

if use_scope
    para0 = [57.9942 -17.4853 -101.6752 -103.2967 -285.1151 20.2564*ones(1, length(ind_mktsize)+1) zeros(1, length(ind_mktsize_scope)+1)];
elseif est_corr
    para0 = [16.2776  -31.2948  -38.1102  -62.1981 -162.3070   14.6300   24.4327   94.0368 0];
else
    para0 = [16.2776  -31.2948  -38.1102  -62.1981 -162.3070   14.6300   24.4327   94.0368];
end


clear fc_iv0 fc_iv0_tr fc_x prod_attributes endo_dev_ind

objga_fit = @(para) as10_market(para, entry_sim_tr, fc_x_tr, ...
    dev_profit_all_tr, dev_profit_one_tr,  ...
    ind_para_fc_x, ind_para_sig, ind_para_scope, ind_para_corr, ind_mktsize, N_dev, g_ar_tr,...
    [], [], [], ...
    [], [], ...
    use_scope, brewer_entry_tr, max_profit_top_vec, min_profit_top_vec, ...
    ind_mktsize_scope, [], [], ...
    brewer_fc_x_tr, brewer_g_ar_tr, brewer_g_ar_tr_boot, N_brewer_dev, fc_err, ...
    est_corr, entry_dev2, N_dev2, corr_sim, dev2_ind_f1_tr, dev2_ind_f2_tr, corr_g_ar, [],...
    reps, sig_level);


objga = @(para) as10_market(para, entry_sim_tr, fc_x_tr, ...
    dev_profit_all_tr, dev_profit_one_tr,  ...
    ind_para_fc_x, ind_para_sig, ind_para_scope, ind_para_corr, ind_mktsize, N_dev, g_ar_tr,...
    entry_sim_boot, fc_x_boot, dev_profit_all_boot, ...
    dev_profit_one_boot, g_ar_boot, ...
    use_scope, brewer_entry_tr, max_profit_top_vec, min_profit_top_vec, ...
    ind_mktsize_scope, max_profit_top_vec_boot, min_profit_top_vec_boot, ...
    brewer_fc_x_tr, brewer_g_ar_tr, brewer_g_ar_tr_boot, N_brewer_dev, fc_err, ...
    est_corr, entry_dev2, N_dev2, corr_sim, dev2_ind_f1_tr, dev2_ind_f2_tr, corr_g_ar, corr_g_ar_boot,...
    reps, sig_level);

%% step 1. find the minima of the test statistic
para_test = nan(N_test, length(para0));
fval_test = nan(N_test, 1);
lb = -2000*ones(1, length(para0));
ub =  2000*ones(1, length(para0));

if est_corr
    lb(ind_para_corr) = 0;
end

parfor k = 1 : N_test    
    if ~ispc
        options = optimoptions('surrogateopt','InitialPoints',para0 + randn/10, 'UseParallel', true, 'MaxFunctionEvaluations',3000);
    else
        options = optimoptions('surrogateopt','InitialPoints',para0 + randn/10, 'UseParallel', false, 'MaxFunctionEvaluations',3000);
    end
    [para_test(k, :), fval_test(k)] = surrogateopt(objga_fit, lb, ub, options);
end

[~, mind] = min(fval_test);
para_pa = para_test(mind, :);
[T_para_pa, cv_pa] = objga(para_pa);

save('minimizer', 'para_pa', 'cv_pa', 'T_para_pa', 'use_scope', 'est_corr');
clear objga_fit lb ub options

%% Step 2. consider small disturbances to the minima obtained in Step 1
para_rand0 = randn(N_perturb, length(para_pa), N_grid)/5;
para_gms_cell0 = cell(N_grid, 1);

parfor np = 1 : N_grid
    para_gms_cell0{np} = gen_gms_perturb(para_pa,...
        para_rand0(:, :, np), N_perturb, objga);
end

para_gms0 = [];
for np = 1 : N_grid
    para_gms0 = [para_gms0; para_gms_cell0{np}];
end
clear para_rand0 para_gms_cell0

save('gms_cs0', 'para_gms0', 'use_scope', 'est_corr')

%% step 3. consider big disturbances to a set of para values obtained in Step 2
if ~isempty(para_gms0)
    % sample a set of para values obtained in step 2
    para_bd = [];
    for k = 1 : length(para_pa)
        [~, min_ind_k] = min(para_gms0(:, k));
        [~, max_ind_k] = max(para_gms0(:, k));
        para_bd = [para_bd; para_gms0(min_ind_k, :)];
        para_bd = [para_bd; para_gms0(max_ind_k, :)];
    end


    para_gms_select = [para_bd(randsample(size(para_bd, 1), N_grid, 'true'), :);...
        para_gms0(randsample(size(para_gms0, 1), N_grid, 'true'), :)];
    clear para_bd


    % consider big disturbances to the sampled set of para values obtained in Step 2
    para_rand = randn(N_perturb, length(para_pa), N_grid*2);

    parfor np = 1 : N_grid*2
        para_gms_cell{np} = gen_gms_perturb(para_gms_select(np, :),...
            para_rand(:, :, np), N_perturb, objga);
    end
    clear para_rand


    para_gms = [];
    for np = 1 : N_grid*2
        para_gms = [para_gms; para_gms_cell{np}];
    end
    clear para_gms_cell

    % save results
    save([output_folder, 'gms_cs'], 'para_gms', 'para_pa', 'use_scope', 'est_corr');
end

